
*
clear
import delimited "C:\...\PES_data_anonymised.csv", bindquote(strict) case(preserve) 
destring, replace
save PES_data_anonymised, replace

** Figure 1: Top motivations to maintain shade trees on farms
clear
use PES_data_anonymised

gr bar (mean) agrf_mot_shade_cocoa agrf_mot_productivity agrf_mot_soil_moisture agrf_mot_timber agrf_mot_building agrf_mot_soil_fertility agrf_mot_firewood agrf_mot_medicine agrf_mot_shade, legend(label(1 "cocoa shade") label(2 "productivity") label(3 "soil moisture") label(4 "timber") label(5 "building material") label(6 "soil fertility") label(7 "firewood") label(8 "medicine") label(9 "work in shade") )  graphregion(color(white))



** Figure 2: Violin plots of outcomes

clear
use PES_data_anonymised

gen t_group=.
replace t_group=1 if vil_treated==0 & SC==0
replace t_group=2 if vil_treated==0 & SC==1
replace t_group=3 if vil_treated==1 & SC==0
replace t_group=4 if vil_treated==1 & SC==1
label define gr 1 "NT nSC" 2 "NT SC" 3 "T nSC" 4 "T SC"
label values  t_group gr


gen t_group2=.
replace t_group2=1 if vil_treated==0 & SC==0
replace t_group2=2 if vil_treated==0 & SC==1
replace t_group2=3 if vil_treated==1 & SC==0
replace t_group2=4 if vil_treated==1 & SC==1 & programme=="Carbone 2122"
replace t_group2=5 if vil_treated==1 & SC==1  & programme=="Carbone 2223"
label define gr1 1 "NT nSC" 2 "NT SC" 3 "T nSC" 4 "T SC" 5 "T2 SC"
label values  t_group2 gr1

gen t_group3=.
replace t_group3=1 if vil_treated==0 & SC==0
replace t_group3=2 if vil_treated==0 & SC==1
replace t_group3=3 if vil_treated==1 & SC==0
replace t_group3=4 if vil_treated==1 & SC==1 & Yeartype=="21/22"
replace t_group3=5 if vil_treated==1 & SC==1  & Yeartype=="21/22 & 22/23"
replace t_group3=6 if vil_treated==1 & SC==1  & Yeartype=="22/23"
label define gr2 1 "NT nSC" 2 "NT SC" 3 "T nSC" 4 "T1 SC" 5 "T1&2 SC" 6 "T2 SC"
label values  t_group3 gr2


violinplot density_ac1, over(t_group) vertical graphregion(color(white)) title("Naturally occurring") ytitle("")
gr save gr0.gph, replace
violinplot density_ac2, over(t_group) vertical graphregion(color(white)) title("Seedlings") ytitle("")
gr save gr1.gph, replace

violinplot density_ac3, over(t_group) vertical graphregion(color(white)) title("Wildlings") ytitle("")
gr save gr2.gph, replace

violinplot density_tot, over(t_group) vertical graphregion(color(white)) title("Total") ytitle("")
gr save gr3.gph, replace

violinplot density_kill if density_kill!=0, over(t_group) vertical graphregion(color(white)) title("Removal") ytitle("")
gr save gr4.gph, replace


graph combine gr0.gph gr1.gph gr3.gph gr4.gph  , col(2) graphregion(color(white))



** Table SI 1
clear

use PES_data_anonymised
collapse (mean) shade_mean cocoa_area urban_area_share def10_19 forest_cover2019 vil_treated , by(iv)


local vars shade_mean cocoa_area urban_area_share def10_19 forest_cover2019 



{
* Initialize row counter for Excel export
local row = 2

* Create an Excel file with headers
putexcel set summary_statistics_village_t.xlsx, replace
putexcel A1=("Variable") B1=("Grand Mean") C1=("N") ///
J1=("Group 1 Mean") K1=("Group 2 Mean") L1=("Difference") M1=("P-value") N1=("N1") O1=("N2")  , bold 

* Loop through each variable
foreach var of local vars {

 local var_label : variable label `var'

    * Calculate summary statistics for the whole sample (optional)
    summarize `var'
  local grandmean = r(mean) 
  local gransample = r(N)
  local grandsd= r(sd)

  

	
	{
	    * Calculate mean differences and t-test p-values - vil_treated
    ttest `var', by(vil_treated)
    * Store the results in macros for exporting
    local dmean1 = r(mu_1)
    local dmean2 = r(mu_2)
    local ddiff = `dmean2' - `dmean1'
    local dpvalue = r(p)
	local dN_1 = r(N_1)
    local dN_2 = r(N_2)
	
	
	local dsd1 = r(sd_1)
	local dsd2 = r(sd_2)
	local ddiff_se = r(se)
	}
	

	
	

    * Export the results to Excel
    putexcel A`row'="`var'" B`row'=`grandmean' C`row'=`gransample' ///
	J`row'=`dmean1' K`row'=`dmean2' L`row'=`ddiff' M`row'=`dpvalue' N`row'=`dN_1' O`row'=`dN_2' ///
	
    
    * Increment row counter
    local row = `row' + 1
	
	* standard errors
	putexcel B`row'=`grandsd' ///
	J`row'=`dsd1' K`row'=`dsd2' L`row'=`ddiff_se' ///
	
    
	* Increment row counter
    local row = `row' + 1
	
	
}

* Save the Excel file
putexcel save 

}









** Table SI2  Household and plot level summary statistics

clear
use PES_data_anonymised


gen t_group=.
replace t_group=1 if vil_treated==0 & SC==0
replace t_group=2 if vil_treated==0 & SC==1
replace t_group=3 if vil_treated==1 & SC==0
replace t_group=4 if vil_treated==1 & SC==1


local vars cocoaplots_dist_vc_pp_KM pde_established_year_pp  plu_primary_forest  cv_Ghana cv_Francais cv_Bresil    cocoaplots_map_pp acp_documentation_yn_pp owner density_ac1 density_ac2 density_ac3  density_tot2  density_kill tot_area ppi hhh_age      member  leader    knowledge training   crd_constraint  buyer_premium live_onfarm wf_manager_n n_input farm_challenge_pest_c yield1 


{
* Initialize row counter for Excel export
local row = 2

* Create an Excel file with headers
putexcel set summary_statistics_plot2_t.xlsx, replace
putexcel A1=("Variable") B1=("Grand Mean") C1=("N") ///
D1=("Group 1 Mean") E1=("Group 2 Mean") F1=("Difference") G1=("P-value") H1=("N1") I1=("N2") ///
J1=("Group 1 Mean") K1=("Group 2 Mean") L1=("Difference") M1=("P-value") N1=("N1") O1=("N2") ///
P1=("Group 1 Mean") Q1=("Group 2 Mean") R1=("Difference") S1=("P-value") T1=("N1") U1=("N2") ///        , bold 

* Loop through each variable
foreach var of local vars {

 local var_label : variable label `var'

    * Calculate summary statistics for the whole sample (optional)
    summarize `var'
  local grandmean = r(mean) 
  local gransample = r(N)
  local grandsd= r(sd)

  
  {
    * Calculate mean differences and t-test p-values - SC
    ttest `var', by(SC)
    * Store the results in macros for exporting
    local mean1 = r(mu_1)
    local mean2 = r(mu_2)
    local diff = `mean2' - `mean1'
    local pvalue = r(p)
	local N_1 = r(N_1)
    local N_2 = r(N_2)

	
	local sd1 = r(sd_1)
	local sd2 = r(sd_2)
	local diff_se = r(se)

	}
	
	{
	    * Calculate mean differences and t-test p-values - vil_treated
    ttest `var', by(vil_treated)
    * Store the results in macros for exporting
    local dmean1 = r(mu_1)
    local dmean2 = r(mu_2)
    local ddiff = `dmean2' - `dmean1'
    local dpvalue = r(p)
	local dN_1 = r(N_1)
    local dN_2 = r(N_2)
	
	
	local dsd1 = r(sd_1)
	local dsd2 = r(sd_2)
	local ddiff_se = r(se)
	}
	
		{
	    * Calculate mean differences and t-test p-values - treated INDIVIDUAL
    ttest `var' if SC==1 & vil_treated==1, by(treated)
    * Store the results in macros for exporting
    local tmean1 = r(mu_1)
    local tmean2 = r(mu_2)
    local tdiff = `tmean2' - `tmean1'
    local tpvalue = r(p)
	local tN_1 = r(N_1)
    local tN_2 = r(N_2)

	
	local tsd1 = r(sd_1)
	local tsd2 = r(sd_2)
	local tdiff_se = r(se)
	}
	
	
	
	

    * Export the results to Excel
    putexcel A`row'="`var'" B`row'=`grandmean' C`row'=`gransample' ///
	D`row'=`mean1' E`row'=`mean2' F`row'=`diff' G`row'=`pvalue' H`row'=`N_1' I`row'=`N_2' ///
	J`row'=`dmean1' K`row'=`dmean2' L`row'=`ddiff' M`row'=`dpvalue' N`row'=`dN_1' O`row'=`dN_2' ///
	P`row'=`tmean1' Q`row'=`tmean2' R`row'=`tdiff' S`row'=`tpvalue' T`row'=`tN_1' U`row'=`tN_2' ///
    
    * Increment row counter
    local row = `row' + 1
	
	* standard errors
	putexcel B`row'=`grandsd' D`row'=`sd1' E`row'=`sd2' F`row'=`diff_se' ///
	J`row'=`dsd1' K`row'=`dsd2' L`row'=`ddiff_se' ///
	P`row'=`tsd1' Q`row'=`tsd2' R`row'=`tdiff_se' ///
    
	* Increment row counter
    local row = `row' + 1
	
	
}

* Save the Excel file
putexcel save 

}




*** table SI 3: Summay statistics about PBC, motivations and satisfaction



clear
use PES_data_anonymised


foreach var in rct_value_2 rct_value_4 rct_value_6 rct_value_10 rct_value_8 rct_value_1_ rct_value_3 rct_value_5 rct_value_7 rct_value_9 rct_value_11 rct_value_12 rct_value_13 {
replace `var'= 6- `var'
}


local vars   rct_value_2 rct_value_4 rct_value_6 rct_value_10 rct_value_8 rct_value_1_ soc_norm_agrof  agrofor_satisfaction agrofor_benefits  fut_farm_remove fut_farm_less fut_farm_more fut_farm_change  fut_farm_pruning fut_farm_same fut_farm_more_fruit  fut_farm_more_timber  rct_value_3 rct_value_5 rct_value_7 rct_value_9 rct_value_11 rct_value_12 rct_value_13 agrf_mot_shade_cocoa agrf_mot_productivity	agrf_mot_shade	agrf_mot_buyers_require	agrf_mot_supply	agrf_mot_timber	agrf_mot_soil_moisture	agrf_mot_medicine	agrf_mot_firewood	agrf_mot_building	agrf_mot_soil_fertility	agrf_mot_labour	agrf_mot_existence	agrf_mot_everyone_has	agrf_mot_production_cost	agrf_mot_project	agrf_mot_income	agrf_mot_recognition	agrf_mot_revenue	agrf_mot_fruit	agrf_mot_fodder	agrf_mot_stake_yam	agrf_mot_no_purpose	agrf_mot_other



{
* Initialize row counter for Excel export
local row = 2

* Create an Excel file with headers
putexcel set summary_statistics_village_likert_t.xlsx, replace
putexcel A1=("Variable") B1=("Grand Mean") C1=("N") ///
D1=("Group 1 Mean") E1=("Group 2 Mean") F1=("Difference") G1=("P-value") H1=("N1") I1=("N2") ///
J1=("Group 1 Mean") K1=("Group 2 Mean") L1=("Difference") M1=("P-value") N1=("N1") O1=("N2") ///
P1=("Group 1 Mean") Q1=("Group 2 Mean") R1=("Difference") S1=("P-value") T1=("N1") U1=("N2") ///        , bold 

* Loop through each variable
foreach var of local vars {

 local var_label : variable label `var'

    * Calculate summary statistics for the whole sample (optional)
    summarize `var'
  local grandmean = r(mean) 
  local gransample = r(N)
  local grandsd= r(sd)

  
  {
    * Calculate mean differences and t-test p-values - SC
    ttest `var', by(SC)
    * Store the results in macros for exporting
    local mean1 = r(mu_1)
    local mean2 = r(mu_2)
    local diff = `mean2' - `mean1'
    local pvalue = r(p)
	local N_1 = r(N_1)
    local N_2 = r(N_2)

	
	local sd1 = r(sd_1)
	local sd2 = r(sd_2)
	local diff_se = r(se)

	}
	
	{
	    * Calculate mean differences and t-test p-values - vil_treated
    ttest `var', by(vil_treated)
    * Store the results in macros for exporting
    local dmean1 = r(mu_1)
    local dmean2 = r(mu_2)
    local ddiff = `dmean2' - `dmean1'
    local dpvalue = r(p)
	local dN_1 = r(N_1)
    local dN_2 = r(N_2)
	
	
	local dsd1 = r(sd_1)
	local dsd2 = r(sd_2)
	local ddiff_se = r(se)
	}
	
		{
	    * Calculate mean differences and t-test p-values - treated INDIVIDUAL
    ttest `var' if SC==1 & vil_treated==1, by(treated)
    * Store the results in macros for exporting
    local tmean1 = r(mu_1)
    local tmean2 = r(mu_2)
    local tdiff = `tmean2' - `tmean1'
    local tpvalue = r(p)
	local tN_1 = r(N_1)
    local tN_2 = r(N_2)

	
	local tsd1 = r(sd_1)
	local tsd2 = r(sd_2)
	local tdiff_se = r(se)
	}
	
	
	
	

    * Export the results to Excel
    putexcel A`row'="`var'" B`row'=`grandmean' C`row'=`gransample' ///
	D`row'=`mean1' E`row'=`mean2' F`row'=`diff' G`row'=`pvalue' H`row'=`N_1' I`row'=`N_2' ///
	J`row'=`dmean1' K`row'=`dmean2' L`row'=`ddiff' M`row'=`dpvalue' N`row'=`dN_1' O`row'=`dN_2' ///
	P`row'=`tmean1' Q`row'=`tmean2' R`row'=`tdiff' S`row'=`tpvalue' T`row'=`tN_1' U`row'=`tN_2' ///
    
    * Increment row counter
    local row = `row' + 1
	
	* standard errors
	putexcel B`row'=`grandsd' D`row'=`sd1' E`row'=`sd2' F`row'=`diff_se' ///
	J`row'=`dsd1' K`row'=`dsd2' L`row'=`ddiff_se' ///
	P`row'=`tsd1' Q`row'=`tsd2' R`row'=`tdiff_se' ///
    
	* Increment row counter
    local row = `row' + 1
	
	
}

* Save the Excel file
putexcel save 

}





******************
** Analysis*******
******************


* load data
clear
use PES_data_anonymised

* define covariate sets
global plot_cov "cocoaplots_dist_vc_pp_KM pde_established_year_pp  plu_primary_forest  cv_Ghana cv_Francais cv_Bresil    i.cocoaplots_map_pp i.acp_documentation_yn_pp owner  "
global indhh_cov "  ppi hhh_age    leader    knowledge training   i.crd_constraint  i.buyer_premium live_onfarm wf_manager_n n_input farm_challenge_pest_c   "
global vil_cov "shade_mean cocoa_area urban_area def10_19 forest_cover2019 " 

gen vil_treated_SC=0
replace vil_treated_SC=1 if vil_treated==1 & SC==1

** Table 4 SI:  Summary results on shade tree density, main model
churdle exp density_ac2 1.vil_treated i.SC vil_treated_SC  shade_mean   i.u [pw=weight] ,  select(1.vil_treated i.SC vil_treated_SC   shade_mean    i.u) ll(0) vce(cl iv) difficult iterate(400)
 outreg2 using myregnofewpl.doc, replace ctitle(`var') addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated i.SC vil_treated_SC)
 margins, dydx(vil_treated_SC) post
estimates store m_density_ac2
 
foreach var in  density_ac1   density_kill density_tot2   {
 churdle exp `var' 1.vil_treated i.SC vil_treated_SC   shade_mean  i.u [pw=weight] ,  select(1.vil_treated i.SC vil_treated_SC  shade_mean   i.u) ll(0) vce(cl iv) difficult iterate(400)
outreg2 using myregnofewpl.doc, append ctitle(`var') addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p)) keep(1.vil_treated i.SC vil_treated_SC)
 margins, dydx(vil_treated_SC) post
 estimates store m_`var'
}
 churdle exp density_ac3 1.vil_treated i.SC vil_treated_SC shade_mean  [pw=weight] ,  select(1.vil_treated i.SC vil_treated_SC shade_mean  ) ll(0) vce(cl iv) difficult iterate(400)
outreg2 using myregnofewpl.doc, append ctitle(`var') addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p)) keep(1.vil_treated i.SC vil_treated_SC )
 margins, dydx(vil_treated_SC) post
 estimates store m_density_ac3
 

 ** Figure 3: Marginal intent to treat effect (ITT) on tree density by tree origin over two years, 95% CI.
 
label var density_ac1 "Native"
label var density_ac2 "Seedlings"
label var density_ac3 "Wildlings"
label var density_kill "Removal"
label var density_tot2 "Overall"
 
 

coefplot m_density_ac2 ///
         m_density_ac3 ///	 
         m_density_ac1 ///
		 m_density_kill ///
		 m_density_tot2, ///
    keep(vil_treated_SC) ///
    xline(0, lpattern(dash)) ///
    swapnames aseq  ///
	bylabels("") ///
    title("") ///
    ylabel(, angle(horizontal)) ///
    ytitle("") ///
    xtitle("ITT (trees per hectare)") ///
    legend(off) ///
	graphregion(color(white)) ///
	coeflabel(m_density_ac1 = "Native" m_density_ac2 = "Seedlings" m_density_ac3 = "Wildlings" m_density_kill = "Removal" m_density_tot2 = "Overall") ///
	ciopts(recast(rcap) color(black))   msymbol(O) mcolor(black)   citop ///
	    mlabel(@b) ///
    mlabposition(2) ///
    mlabcolor(black) ///
   format(%5.2f)
	





** Tables SI 5-9: Test robustness for each agroforestry outcome
foreach var in  density_ac1 density_ac2 density_tot2 density_kill    {

capture noisily: reg `var' 1.vil_treated##i.SC i.u [pw=weight], vce(cl iv)  
capture noisily: outreg2 using myregnofew`var'.doc, replace ctitle(OLS) addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: reg `var' 1.vil_treated##i.SC i.u shade_mean [pw=weight], vce(cl iv) 
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(OLS_meandens) addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: reg `var' 1.vil_treated##i.SC i.u i.iv [pw=weight], vce(cl iv) 
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(OLS_fe) addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC   i.u [pw=weight] ,  select(1.vil_treated##i.SC   i.u) ll(0) vce(cl iv) difficult iterate(400)
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_1) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))   keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC shade_mean  i.u [pw=weight],  select(1.vil_treated##i.SC shade_mean i.u) ll(0) vce(cl iv) difficult iterate(400) // increases the proportion of adopters, but not the n of density
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_dens) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))   keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC i.iv  i.u [pw=weight],  select(1.vil_treated##i.SC i.iv i.u) ll(0) vce(cl iv) difficult iterate(400)
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_fe) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC ${plot_cov} shade_mean i.u [pw=weight],  select(1.vil_treated##i.SC ${plot_cov} shade_mean i.u) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control1) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC ${indhh_cov} shade_mean  i.u [pw=weight],  select(1.vil_treated##i.SC ${indhh_cov} shade_mean i.u) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control2) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

capture noisily: churdle exp `var' 1.vil_treated##i.SC ${vil_cov} shade_mean  i.u [pw=weight],  select(1.vil_treated##i.SC ${vil_cov}  shade_mean i.u) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control3) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

capture noisily: churdle exp `var' 1.vil_treated##i.SC ${plot_cov} ${vil_cov} shade_mean  i.u [pw=weight],  select(1.vil_treated##i.SC  ${plot_cov} ${vil_cov} shade_mean i.u) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control_4) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)


capture noisily: churdle exp `var' 1.vil_treated##i.SC ${indhh_cov} ${plot_cov} ${vil_cov} shade_mean  i.u [pw=weight],  select(1.vil_treated##i.SC ${indhh_cov} ${plot_cov} ${vil_cov} shade_mean i.u) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control_all) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

}


foreach var in  density_ac3    {
capture noisily: reg `var' 1.vil_treated##i.SC i.u [pw=weight], vce(cl iv)  
capture noisily: outreg2 using myregnofew`var'.doc, replace ctitle(OLS) addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: reg `var' 1.vil_treated##i.SC i.u shade_mean [pw=weight], vce(cl iv) 
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(OLS_meandens) addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: reg `var' 1.vil_treated##i.SC i.u i.iv [pw=weight], vce(cl iv) 
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(OLS_fe) addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC    [pw=weight] ,  select(1.vil_treated##i.SC   ) ll(0) vce(cl iv) difficult iterate(400)
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_1) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))   keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC shade_mean   [pw=weight],  select(1.vil_treated##i.SC shade_mean ) ll(0) vce(cl iv) difficult iterate(400) // increases the proportion of adopters, but not the n of density
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_dens) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))   keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC i.iv   [pw=weight],  select(1.vil_treated##i.SC i.iv ) ll(0) vce(cl iv) difficult iterate(400)
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_fe) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC ${plot_cov} shade_mean [pw=weight],  select(1.vil_treated##i.SC ${plot_cov} shade_mean ) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control1) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
capture noisily: churdle exp `var' 1.vil_treated##i.SC ${indhh_cov} shade_mean   [pw=weight],  select(1.vil_treated##i.SC ${indhh_cov} shade_mean ) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control2) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

capture noisily: churdle exp `var' 1.vil_treated##i.SC ${vil_cov} shade_mean   [pw=weight],  select(1.vil_treated##i.SC ${vil_cov}  shade_mean ) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control3) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

capture noisily: churdle exp `var' 1.vil_treated##i.SC ${plot_cov} ${vil_cov} shade_mean  [pw=weight],  select(1.vil_treated##i.SC  ${plot_cov} ${vil_cov} shade_mean ) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control_4) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)


capture noisily: churdle exp `var' 1.vil_treated##i.SC ${indhh_cov} ${plot_cov} ${vil_cov} shade_mean   [pw=weight],  select(1.vil_treated##i.SC ${indhh_cov} ${plot_cov} ${vil_cov} shade_mean ) ll(0) vce(cl iv) difficult iterate(400) // adding plot level controls
capture noisily: outreg2 using myregnofew`var'.doc, append ctitle(DH_control_all) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

}



** Table 10:  Impact on diversity of planted seedlings – OLS regression

foreach var in tree_diversity {

reg `var' 1.vil_treated##i.SC i.u [pw=weight], vce(cl iv)
 outreg2 using myreg`var'.doc, replace ctitle(OLS) addstat(Adjusted R2, e(r2_a)) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
reg `var' 1.vil_treated##i.SC i.u i.iv [pw=weight], vce(cl iv) 
 outreg2 using myreg`var'.doc, append ctitle(OLS_fe) addstat(Adjusted R2, e(r2_a)) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

reg `var' 1.vil_treated##i.SC ${plot_cov} i.u [pw=weight], vce(cl iv) 
 outreg2 using myreg`var'.doc, append ctitle(OLS_control) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

reg `var' 1.vil_treated##i.SC ${indhh_cov} i.u  [pw=weight], vce(cl iv) 
 outreg2 using myreg`var'.doc, append ctitle(OLS_control2) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
 
reg `var' 1.vil_treated##i.SC ${vil_cov} i.u  [pw=weight], vce(cl iv) 
 outreg2 using myreg`var'.doc, append ctitle(OLS_control2) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

reg `var' 1.vil_treated##i.SC ${plot_cov} ${indhh_cov} i.u ${vil_cov} [pw=weight], vce(cl iv) 
 outreg2 using myreg`var'.doc, append ctitle(OLS_control3) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)

}


**Table 4 and Table SI 11: Impacts on the elements of the TPB – OLS regression

* rescaling likert outcomes
foreach var in  rct_value_2 rct_value_4 rct_value_6  rct_value_8  soc_norm_agrof  agrofor_satisfaction agrofor_benefits {
replace `var'= 6- `var'
}
sum rct_value_2 rct_value_4 rct_value_6  rct_value_8  soc_norm_agrof  agrofor_satisfaction agrofor_benefits

reg rct_value_2 1.vil_treated##i.SC  i.u  [pw=weight], vce(cl iv)  
outreg2 using myregs.doc, replace ctitle(rct_value_2) addstat(Adjusted R2, e(r2_a))
foreach var in   rct_value_4 rct_value_6  rct_value_8  soc_norm_agrof agrofor_satisfaction agrofor_benefits  {
reg `var' 1.vil_treated##i.SC i.u  [pw=weight], vce(cl iv)  
outreg2 using myregs.doc, append ctitle(`var') addstat(Adjusted R2, e(r2_a))
}


** Table SI 12: Robustness check of the impact on the elements of the TPB – OLS regression
reg rct_value_2 1.vil_treated##i.SC ${plot_cov} ${indhh_cov} i.u  [pw=weight], vce(cl iv)  
outreg2 using myregs.doc, replace ctitle(rct_value_2) addstat(Adjusted R2, e(r2_a))
foreach var in   rct_value_4 rct_value_6  rct_value_8  soc_norm_agrof agrofor_satisfaction agrofor_benefits  {
reg `var' 1.vil_treated##i.SC ${plot_cov} ${indhh_cov} i.u  [pw=weight], vce(cl iv)  
outreg2 using myregs.doc, append ctitle(`var') addstat(Adjusted R2, e(r2_a))
}


** Table SI 13:  Impacts on tree planting motivations – OLS regression 
egen nmot=rowtotal(agrf_mot_*)
reg nmot  1.vil_treated##i.SC   i.u [pw=weight], r
 outreg2 using myregmotivation2.doc, replace ctitle(nmot) addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p)) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
foreach var in shade_cocoa	productivity	shade	buyers_require	supply	timber	soil_moisture	medicine	firewood	building	soil_fertility	labour	existence	everyone_has	production_cost	project	income	recognition	revenue	fruit	fodder	stake_yam	no_purpose	other{

reg agrf_mot_`var'  1.vil_treated##i.SC   i.u [pw=weight], r
 outreg2 using myregmotivation2.doc, append ctitle(`var') addstat(Adjusted R2, e(r2_a), Chi2 p value, e(p)) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
}



** Table SI 14: Impacts on future vision - OLS regression
reg fut_farm_remove 1.vil_treated##i.SC i.u [pw=weight], r 
capture: outreg2 using myregfut.doc, replace ctitle(fut_farm_remove) addstat(Adjusted R2, e(r2_a), Chi 2 p value, e(p)) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
foreach var in  fut_farm_less fut_farm_more fut_farm_change  fut_farm_pruning fut_farm_same fut_farm_more_fruit  fut_farm_more_timber {

reg `var' 1.vil_treated##i.SC i.u [pw=weight], r  
outreg2 using myregfut.doc, append ctitle(`var') addstat(Adjusted R2, e(r2_a), Chi 2 p value, e(p)) keep(1.vil_treated 1.SC 1.vil_treated#1.SC)
}



** Table 5 and Table SI 15: Mediation analysis of PES through the elements of TPB
churdle exp density_ac2 rct_value_2 rct_value_4 rct_value_6  rct_value_8  shade_mean   i.u [pw=weight] ,  select( rct_value_2 rct_value_4 rct_value_6  rct_value_8  shade_mean    i.u) ll(0) vce(cl iv) difficult iterate(400)
 outreg2 using myregnofewpl.doc, replace ctitle(TPB) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p))  
 
churdle exp density_ac2 soc_norm_agrof shade_mean  i.u [pw=weight] ,  select(soc_norm_agrof  shade_mean   i.u) ll(0) vce(cl iv) difficult iterate(400)
outreg2 using myregnofewpl.doc, append ctitle(TPB-PCA) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p)) 
 
 churdle exp density_ac2 1.vil_treated##i.SC   shade_mean  i.u [pw=weight] ,  select(1.vil_treated##i.SC  shade_mean   i.u) ll(0) vce(cl iv) difficult iterate(400)
outreg2 using myregnofewpl.doc, append ctitle(PES) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p)) 

 churdle exp density_ac2 1.vil_treated##i.SC  rct_value_2 rct_value_4 rct_value_6  rct_value_8 shade_mean  i.u [pw=weight] ,  select(1.vil_treated##i.SC rct_value_2 rct_value_4 rct_value_6  rct_value_8 shade_mean   i.u) ll(0) vce(cl iv) difficult iterate(400)
outreg2 using myregnofewpl.doc, append ctitle(PES - TPB) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p)) 

 churdle exp density_ac2 1.vil_treated##i.SC  soc_norm_agrof shade_mean  i.u [pw=weight] ,  select(1.vil_treated##i.SC  shade_mean soc_norm_agrof  i.u) ll(0) vce(cl iv) difficult iterate(400)
outreg2 using myregnofewpl.doc, append ctitle(PES - TPB PCA) addstat(Adjusted R2, e(r2_p), Chi2 p value, e(p)) 


